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This paper presents numerical resnlts for the two-dimensional isotropic Kuramoto-Sivashinsky 
equation (KSE) with an additional nonlinear term and a single independent parameter. Surfaces 
generated by this equation exhibit a certain dependence of the average saturated roughness on the 
system size that indicates power-law shape of the surface spectrum for small wave numbers. This 
leads to a conclusion that although cellular surface patterns of dehnite scale dominate in the range 
of short distances, there are also scale-free long-range height variations present in the large systems. 

The dependence of the spectral exponent on the equation parameter gives some Insight into the 
scaling behavior for large systems. 


I. INTRODUCTION 

The Kuramoto-Sivashinsky equation (KSE) in its dimensionless form for some held h can be written as [IHl] 

dth = -V'^h-V'^h+{Vh)'^. ( 1 ) 

This equation stands as a paradigmatic model for chaotic spatially extended systems and can be used to study the 
connections between chaotic dynamics at small scales and apparent stochastic behavior at large scales. It is an example 
of an extended, deterministic dynamical system that exhibits complex spatio-temporal phenomena. The KSE has 
been derived for the purpose of describing the intrinsic instabilities in laminar hame fronts ^ and phase-dynamics 
in reaction-diffusion systems [2]. The equation Q in one- and two-dimensional cases has been a subject of active 
research for about three decades, and its scaling properties have even been an object of some controversy 13 El- 

This paper presents some results obtained from a less researched generalized version of KSE. There have been many 
different generalizations of the KSE used for different purposes. Some of them involve adding damping terms to Q 
[11 IS] , some introduce spatial anisotropy |6] or additive random noise |7] . 

In the case presented here, there is an additional nonlinear term V^(V/i)^ introduced to the KSE Q. The two- 
dimensional generalized Kuramoto-Sivashinsky equation of this form with an additive Gaussian white noise has been 
used as a model equation for amorphous solid surface growth BM- The equation in this model originally has 
five parameters that are needed in order to reproduce the experimental data in the simulations and to examine the 
correspondence with microscopic properties of the surface growth process [S]. However, for theoretical investigations of 
the long-time and large-scale behavior of the system, the noise term can be neglected and the remaining deterministic 
equation can be rescaled into the dimensionless form with only one independent parameter a: 

dth = -V^h - + {Vhf . (2) 

Equation ([^ has also been used as a model for nano-scale pattern formation induced by ion beam sputtering [ini E] ■ 

The two-dimensional generalized Kuramoto-Sivashinsky equation (GKSE) (§ in the context of this paper describes 
the evolution of a (2-|-l)-dimensional surface, i.e., a surface which is defined as a function on a two-dimensional plane 
and is growing in the direction perpendicular to that plane. The surface profile h{r, t) is defined as the surface height 
h at the position r = {x, y) on the square [0, of size L in the plane at time t or, more generally, as a function: 

h : [0, L]2 X M+ ^ K . (3) 

We solve the equation (§ for different values of the parameter a using the finite difference method with periodic 
boundary conditions, the time step At = 0.005, spatial discretization step Aa; = 0.71086127010534. Such a seemingly 
bizarre number for the discretization step Aa; is actually a good approximation of the value that is needed in order 
for the system with periodic boundary conditions to be able to contain ordered patterns that appear in some other 
versions of the generalized KSE. The equation is solved for system sizes L ranging from about 45 to about 1000 (i.e., 
on the N X N lattices with N from 63 to 1400, where L = N Ax) and in two cases up to about 1422 {N = 2000). 
The methods of numerical solution for ([^ are presented and compared in [12]. 


II. KINETICS OF THE SURFACE ROUGHNESS 

The surface roughness wft), also called the surface width, is one of the most important characteristics of a surface 
m- It is defined as the standard deviation, or, synonymously, root mean square (rms) deviation, of the surface height 
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FIG. 1. Kinetics of the surface roughness 'w{t) for short times (t < 500) for surfaces with initial roughness w(0) = 10 evolving 
according to it with parameter values a = 0 {black, dashed line), a — 0.5 {red, continuous line), a = 1 {blue, dotted line). 
System size L « 355 (in lattice units, N = 500). 


h{r,t) from its average value h{t) := {h{r,t))r, at some time t: 


w{t) := ^(^{h{r,t)-h{t)f'^ (4) 

Here and throughout the whole paper we denote the averaging by (• ■ ■) with some subscript that shows over what 
entities the averaging takes place. (•■■)»- denotes spatial average over the whole surface, (• • • )t temporal average, 
{■ ■ •)n ensemble average, (• • • )| 7 -|=r spatial average over all r whose length is r etc. 

For smaller values of parameter a (say, a < 5) the kinetics of w{t) due to the evolution of the surface resulting from 
i]) seems to follow a distinct pattern. Starting from a random surface with some small initial roughness w{t = 0) <C 1, 
the roughness begins to grow almost exponentially, but at some time tx ~ 100 this growth slows down and later on 
crosses over to a stationary regime where it oscillates about its average {saturation) value Wsat (see Fig. [^. 

In the saturation regime, the dynamics of w{t) becomes a statistically stationary process with time independent 
average and other statistical characteristics that are the same for different realizations (different realizations differ in 
the initial surface profile, as the evolution equation itself is deterministic). However, in order to consider the long time 
behavior (after saturation) of a statistically stationary process, it is important to take time averages over sufficiently 
long time intervals, having in mind that the minimal required averaging time must be at least several times longer 
than the correlation time of the process, defined as the time lag value at which the normalized autocorrelation function 
of 'w{t) effectively falls to zero. This correlation time might strongly depend on the parameter a and the system size 
L. 

In the stationary regime the roughness is chaotically oscillating about some average value which we denote Wsat 
and call saturated roughness. This value is calculated as the time average of w{t) in the stationary regime and is 
virtually the same when averaged over different time intervals of the stationary regime (given that these intervals are 
sufficiently long) and for different realizations. 

This saturated surface roughness can be theoretically defined as a time average over an time interval of length T 
that goes to infinity, starting from the time tg where the saturation regime is surely reached: 

Our numerical investigation shows that, for system sizes N and parameter values a considered in this paper, it is 
sufficient to take tg > 3000 and T > 10000 to get the saturation values rcsat that differ less than 5% for different 
realizations. For most of the results presented here, we have used the values iCsat obtained using tg > 9000 and 
T = 11000 and, furthermore, we used the averaged values of several (from 3 to 10) realizations: 

«;:= ((u;(t))^g[9.io3 2.104])^ • (6) 

In this paper we are interested in the surface patterns produced by GKSE ([^ and the dependence of the saturated 
surface roughness on parameter a and system size L. 
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FIG. 2. Left panel: Surfaces h{r, t) (values of h coded in gray-scale) evolving according to § at system size N = 250 {L ~ 178) 
with parameters a = 0, 0.25, 0.5, 1, 2 at time t = 5000 (in the stationary regime). Right panel: Normalized autocorrelation 
functions C{r) as defined in Q of the surfaces that are shown on the left panel. 


III. RESULTING SURFACES 


The surfaces, generated by Eq. ([^ in the stationary regime have a disordered cellular structure with cells whose 
sizes are in a quite narrow interval (see Figs.[2jQ. The resulting patterns seem to have similar appearance for different 
system sizes L (for sizes that are at least several times larger than the typical cell diameter). 

The usual way to investigate the surface patterns is by calculating the surface height correlation function C'(r), 
which is the two-dimensional autocorrelation function of the surface height h{r): 

C{r) = {{h{r') - h) {h{r' + r) - h))^, . (7) 

Since the equation ([^ is isotropic, the correlation functions Q C{r) of the resulting surfaces should statistically 
be independent of the direction of r and, thus, depend only on its absolute value r = |r|. We therefore define the 
isotropic height correlation function C'(r) as averaged over all directions of r: 

C{r) =/{{h{r')-h){h{r'+ r)-h)),\ . (8) 

\ / |r|=r 


Fig. i shows resulting surface patterns for system size L « 178 (or, in lattice units, N = 250) and different values 
of parameter a in and the corresponding normalized height correlation functions §. We see that for the surfaces 
generated by (§ with (2 = 0, corresponding to the KSE 0 case, the height correlation function (|^ has no maximum 
(right panel of Fig. [^, just an area of slower decay at distances r corresponding to the approximate sizes of cells in 
the pattern. For a > 0, the height correlation function obtains a maximum whose distance corresponds to the average 
size of cells. We see that this the distance of this maximum increases with o:, meaning that the cell size increases as 
the parameter a is increased. 

Another thing that can be seen in Fig. [^is that the normalized correlation functions Cir^jw^ decrease slowly for 
small values of a > 0 and faster for larger values. That is the first indication of the influence of parameter ol on 
long-range height correlations. 

Large-scale height variations in surfaces produced by Eq. ([^ become more distinct as the system size is chosen to 
be many times larger than the typical cell size. The examples for some values of a can be seen in Fig. and Fig.[^ 

By repeating the simulations with different system sizes L, one notices a slight dependence of the average saturation 
value of the surface roughness Wsat on the system size. Indeed, the results show that the resulting roughness increases 
when the system size is increased. This shows that the resulting surface profile contains spatial Fourier components 
of ever smaller wave number k (and correspondingly larger wavelength A). Although for small scales (large wave 
numbers) the structures of definite size occur, for larger distances (small wave numbers) we get the height variations 
with long-range dependence, and (as will be shown in the next two sections) this dependence has a scale free character. 
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FIG. 3. Surfaces (values of h coded in gray-scale) at system size N = 1200 {L ~ 853) evolving according to § with parameters 
a = 0 (left panel) and a = 0.25 (right panel) at time t = 6000 (in the stationary regime). 



FIG. 4. Surfaces (values of h coded in gray-scale) at system size N — 1200 {L « 853) evolving according to § with parameters 
a = 0.5 (left panel) and a = 1 (right panel) at time t = 6000 (in the stationary regime). 


IV. POWER-LAW SURFACE SPECTRA AND SCALING OF ROUGHNESS 

This section presents some general considerations about the spatial power-spectral densities (PSD) of surfaces, their 
connection to the surface roughness w, and the effects introduced by the finite system size. The scaling behavior of w 
when the PSD has a power-law shape is derived. These theoretical results are compared to the numerically calculated 
scaling properties of w for the surfaces generated by in the next section. 


A. Surface PSD and roughness 

A two-dimensional surface h(r) is a single valued function on the plane r = {x,y) S In order to avoid the zero 
frequency component in the spectrum, we calculate the spectrum of the surface prohle hc{r) with zero mean: 

hc{r) := h{r) - h, (9) 

where h = {h{r))r is the average height of the surface. Fourier transformation H{k) of this ’centered’ surface profile 
hc{r)-. 


H{k) = 


d^r hc{r) exp(—ifc • r). 


( 10 ) 


Here k = {kx, ky)^ is the wave vector of spatial Fourier components of the surface profile. 

The surface power-spectral density (PSD) is then defined as: 

su = jj^mk)w ( 11 ) 

where L is the size of the segment of the surface analyzed, i.e., r G [0,L]^. The integral of the power spectrum Sk 
over all k is equal to the variance of the surface height cr^ which by our definition is is equal to the square of the 
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surface roughness w^: 


1 


d^k Sk = . 


( 12 ) 


Since our model is isotropic, Sk must depend only on the absolute value |fc|, and we can get the one-dimensional 
power spectrum S{k) of the surface by integrating the two-dimensional power spectrum over all wave vectors k of the 
same absolute value k. This can be done by expressing the wave vector k in the polar coordinates k = {k,(f>), and 
then integrating Sk over all angles (j): 


From this we get that the one-dimensional PSD of the surface can be expressed as: 


S{k) = A (\Hik)\%^^^ = A \Hik)\^ = kSk. 


L2 


(14) 


In Eq. (14) we have used the fact that for isotropic surfaces {\H{k)\'^)\k\=k = \H{k)\'^. The integral of this one¬ 
dimensional PSD S{k) over all wave numbers k also equals to the square of surface roughness: 


1 

27r 


Jdk S{k) = w'^ . 


(15) 


B. Effects of finite system size 

The surfaces in numerical simulations are represented on a (TV x N) matrix that sets limits to the smallest and the 
largest possible wave numbers fcmin and fcmax that can fit into the system, thus, ’filtering’ the theoretically defined 
PSD S{k). If Aa: is the spatial step size in the simulation, then the minimal distinguishable wavelength rmin in the 
system is approximately equal to this discretization step doubled: 

rmin « 2Aa;, (16) 

and maximal wavelength that can fit into the system is of about double system size: 

J'max « 2L = 2NAx . (17) 

Since the distance r corresponds to the wave number fc = ^, we get the minimal and maximal wave numbers for the 
system: 


h ■ 

'^miri 


NAx ’ 


Aa: ' 


(18) 


The square of the numerically calculated surface roughness, expressed according to (15) should then be 


- ikm ■ 


(19) 


If we keep the discretization step Aa: constant, assuming that surface patterns for systems of different sizes (up to 
the smallest wave numbers allowed by the system size) are statistically the same, then, by increasing the system size N, 
according to (18), we reduce the minimal wavenumber fcmin in (19). Therefore the calculated surface roughness must 


grow with the system size and its scaling behavior w{N) when N is increased should be able to give us information 
about the shape of the surface PSD S{k) for small wave numbers fc —)■ 0. 
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C. Power-law spatial spectrum for small wave numbers 


Let us assume that the one-dimensional spatial PSD S{k) ( (l4| ) of a surface has a power-law dependence on k for 
small wave numbers, smaller than some wavelength ks- 


Qfi,\ _ / ^ ^ ^ k < ks, 

^ ' 1 Si{k) for k > ks ■ 


( 20 ) 


Here C is some constant, Si{k) is the shape of the PSD spectrum for high wave numbers k that does not interest 
us, since we are interested in large scale behavior of the system and, furthermore, assume that Si{k) doesn’t change 
when the system size is increased. We also assume that the power for small wave numbers 7 < 1. 

Then, for a system of finite size and fc„iin < ks < kmux, the calculated square of the roughness is expressed by: 


w 


dkS{k) = C / dkk-'* + A. 


( 21 ) 


Since ks and 7 come from the model and femax comes from numerical scheme, the only variable is /cmin which is 
inversely proportional to the system size. Therefore the second term in (21) is just a constant which we denoted by 
A: 


A = 


dkSi{k). 


( 22 ) 


There are two qualitatively distinct cases. One is 7 = 1 which would result in infinite w'^ for a system of infinite 
size and another case is 7 < 1 for which even a system of infinite size would have a finite height variance w^. For 
7 = 1 , the square of the surface roughness w from (21 1 grows linearly with the logarithm of the system size N-. 


Cln ( -^1 +A = C\n(N)+B. 
V«min/ 


Here C and B = A — ln(7Vs) with 


Ns = 


ksNx 


(23) 


(24) 


being constant. The last expression in (23) is more useful, since we don’t know the exact value of Ns- We see that 
when 7 = 1 , the surface roughness w goes to infinity for infinite size system {N -^ 00 ). This means that in this case 
the influence of long-range height variations grows with the system size. 

When 7 < 1 , from (21 1 we get the following scaling relation: 


C 


w 


1-7 


- kl-N) + H = d(i - {N/Ns)-^^-^^) + A 


with constants A, as defined above, and D, which defined as 

C 


D = 


1-7 


(i^) 


1-7 


Again, since we don’t know the value of A^g, we express the scaling of the surface roughness as 

« Cl - C2 

with the constants 

Cl = D A ^ 

and 

C2 = D Nl-^ . 

We see that in this case (for 7 < 1) the surface roughness w has a finite value for a system of infinite size: 

Woo := lim w{N) = \f^\. 

AT—>-oo 


(25) 

(26) 

(27) 

(28) 

(29) 


( 30 ) 
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FIG. 5. (log-log scale) Numerically calculated PSD of the surfaces at time t = 10^ produced by ([^ with parameter a = 0 for 
system sizes (in lattice units) N = 500 (blue triangles), N = 250 (red filled circles) and N = 125 (green squares). The black 
dashed line represents power-law fit with exponent 7 = 1.003. 
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FIG. 6 . (log-log scale) Numerically calculated PSD of the surfaces at time t = 10"* produced by ([^ with parameter a = 1 for 
system sizes (in lattice units) N = 500 (blue triangles), N = 250 (red filled circles) and N = 125 (green squares). The black 
dashed line represents power-law fit with exponent 7 = 0.836. 


V. NUMERICAL RESULTS 

In this section we present our numerical results for surfaces generated by (1^: the one-dimensional power-spectral 
densities (PSD) S{k) ( [l4| ) calculated from the autocorrelation function C'(r) for system sizes up to = 500 and 
the scaling of the roughness w which gives shapes of the S{k), based on the considerations of the previous section. 


A. Surface spectra 


Since the equation ([^ is isotropic, the direction of the wave vector k does not matter in statistical description of 
height variations, and in spectral analysis of the surface patterns, the wave number k = |fc| is sufficient to describe 
the occurring spatial modes. We can therefore analyze one-dimensional surface spectra S{k), defined in ( |14[ ). 

We calculate the one-dimensional surface spectrum S{k) that depends only on the wave number k using the Wiener- 
Khinchin theorem which states that the power-spectral density can be obtained from the Fourier transform of the 
autocorrelation function. Although the height correlation function C'(r) ^ depends only on the absolute value r of 
the shift r, it is nevertheless a two-dimensional autocorrelation function of the surface h{r). Thus, by applying the 
two-dimensional Fourier transform on the height correlation function (|^, we get the two-dimensional PSD (11) from 
which we calculate the one-dimensional PSD S{k) using (14). 














The numerically calculated height correlation function for an isotropic surface h{r) has been defined in 

C{r) = {{hc(r') hc{r' + 


|r|=r 


Here, as before, hdr) = h{r) — h is the surface height with its average value subtracted, the so-called ’centered’ 
surface profile. The two-dimensional PSD Sk is calculated according to ( [IT| ) with ( [To| ): 


S'fc = ^ y d^r' J d^r” hc{r') hc{r'')e 


-r') 


By changing the variable r" r' + r and switching the order of integration, we get 


1 


where, according to Q, 


S'*, = / dVe ‘ / d^r' hc{r')hc{r' + r) 


^ J d^r' /ic(r') /ic(r' -b r) = (/ic(r') hc{r' + r))^, = C{r) 


(31) 


(32) 


(33) 


is the height correlation function of the surface h{r). 

Since the surfaces h{r) generated by (1^ are statistically isotropic, the correlation function (33) effectively depends 
only on r = |r| and can be denoted by C\r). Using this, expressing the two-dimensional integration in (32) over r in 
the polar coordinates r = (r, 4>) and integrating over the angular part, we get 


\kr cos (f) 


Sk = J dr rC{r) J d(j)e 

Substituting this into (13) gives the one-dimensional PSD of the surface: 

S{k) = k2TT J dr r C{r) Jo{kr), 
Here Jo{kr) is the Bessel function of the 1st kind: 

1 

Jo{kr) = ^ J d(/>e 


ikr cos (p 


(34) 


(35) 


(36) 


Figures and show the calculated surface spectra of the surfaces generated by ([^ with parameter values a = 0 
(the KSE case) and a = 1, respectively. In each figure, the spectra for systems of sizes (in lattice units) N = 125, 
N = 250 and N = 500 are shown. Each spectrum is obtained by averaging over 6 different realizations. The figures 
confirm the above assumption that increasing the system size does not change the small scale structures, since the 
spectra at large values of k coincide. The appearance of low wave number modes at larger system sizes and the 
approximate power-law behavior of the spectrum S{k) can also be seen. 

Of course, for larger systems, the numerical calculation of the two-dimensional autocorrelation function 0 and 
surface spectrum (35) directly can take a very long time. Therefore, the possibility to obtain S{k) from scaling of the 
surface roughness is very useful. 


B. Scaling of the surface roughness 


We have investigated the dependence of this saturated surface roughness on the size of the model system for different 
values of the parameter a in 0- We investigated how the square of w estimated by 0 depends on the system size 
N (in lattice units) for sizes N = 63, N = 88, N = 125, N = 175, N = 250, N = 350, N = 500, N = 700, N = 1000, 
N = 1400, N = 2000. The results are shown in Fig. [^for a = 0 and a = 0.25, Fig. [^for a = 0.5 and a = 1. 

For a = 0 (the KSE case), the scaling exponent (as introduced in (20l) 7 = 1 fits the results best. Thus the 
dependence of w'^ on N can be approximated by ( |23[ ) and the results should form a straight line in the log-linear plot, 
and, indeed they do as can be seen in Fig.(filled black circles). The blue line in Fig. [^fitted to data for sizes from 
N = 88 to N = 1400. 
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FIG. 7. (log-linear scale) Black symbols: numerical results for surfaces evolving according to © with a = 0 (filled circles) and 
a = 0.25 (filled squares). Time averaged square of the surface roughness plotted as a function of the system size JV (in 
lattice units). Lines: fits of the results for a = 0 (blue line) and a — 0.25 (red line) by (231 and \21\ respectively, corresponding 
to the scaling surface spectrum ( |20[ ). Resulting fit parameters for a = 0.0: B = —1.0935, C = 1.224; for a = 0.25: 7 = 0.991, 
Cl = 67.770, C 2 = 68.284 . 



FIG. 8 . (log-linear scale) Black symbols: numerical results for surfaces evolving according to © with a = 0.5 (filled triangles) 
and a = 1.0 (filled squares).Time averaged square of the surface roughness plotte d as a function of the system size N (in 
lattice units). Lines: fits of the results for a = 0.5 (blue line) and a = 1.0 (red line) by (271 corresponding to the scaling surface 
spectrum (20l. Resulting fit parameters for a = 0.5: 7 = 0.912, Ci = 4.554, C 2 = 4.768; for a = 1.0: 7 = 0.827, Ci = 1.563, 
C 2 = 1.532 . 


For a = 0.25, the scaling exponent 7 « 0.991 fits the results best and the dependence of w'^ on N can be 
approximated by (27). Fig. shows the calculated results (filled black rectangles). There seems to be a good 
agreement between the data and the fitting curve which is fitted to data for sizes from TV = 125 to = 1400 (red 
line in Fig. [^. 

The results for a = 0.5 and a = 1 give scaling exponents are 7 « 0.912 and 7 « 0.827, respectively, and the 
dependence of w'^ on N can be approximated by (27). Fig. shows the calculated results (black triangles for a = 0.5 
and black rectangles a = 1) together with their respective fits (blue and red lines in Fig. . 


VI. SUMMARY AND DISCUSSION 


The generalized Kuramoto-Sivashinsky equation © produces surfaces with disordered cellular patterns (Figs. 
and 1^ . The size of the average size of a cell depends on equation parameter a and constitutes a definite scale in the 
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surface pattern (see the peaks on the correlation functions in However, in larger systems, the long-range height 
variations of different character become apparent (Figs. and^. We investigate these long-range height variations 
for several values of parameter a by calculating the corresponding scaling relations of the surface roughness. 

The square of the surface roughness w which is by definition equal to the variance of the surface height h(r, t), 




can also be expressed as an integral (15) of the one-dimensional power spectral density S{k), given in 


of the surface over all wave numbers k. Since the finite size L = NAx of the system and the discretization step Ax 
in numerical simulations d efin e the approximate lower and upper cut-off values (181 for the possible wave numbers, 
the estimated value of (19) depends on the system size. As indicated by the surface spectra in Figs. and for 
systems of increasing size, the small-scale patterns (corresponding to higher wave numbers k) remain statistically the 
same, and, additionally, new lower wave number modes arise in larger s yste ms. Therefore by increasing the size of the 
system, from the corresponding change in the calculated value of (19), the shape of the spatial power spectrum 
S{k) can be extracted. 

K the behavior of S{k) for small wave numbers k (say when k < kg, for some kg) follows the inverse power-law 
(20), S{k) oc k~'^, with 7 > 0, this indicates that long-range, scale-free height variations are present. There can be 
qualitatively distinct cases for different values of the exponent 7 . If 7 = 1, then we get the scaling relation (23) which 
implies that w will grow indefinitely as the system size goes to infinity. If, on the other hand, 7 < 1, then we get the 
scaling relation (27), therefore w will approach finite value (30) as the system size goes to infinity. 

For the surfaces generated by numerical simulations of the generalized isotropic Kuramoto-Sivashinsky equation 
(i)- we indeed see the indications that the spatial power spectral density S{k) follows a power-law ( [^ , since tte 


theoretically calculated scaling relations (23) and ([^ fit the numerically established scaling of weir~(see Figs. 
and[^. For the Kuramoto-Sivashinsky equation ([^, we get the spectral exponent 7 = 1 and the scaling relation ([^3). 
For parameter values a > 0 the scaling exponent 7 decreases resulting in the scaling relation ( |27| giving the surfaces of 
finite roughness as the size of the system goes to infinity, and thus showing that the generalized Kuramoto-Sivashinsky 
equation ([^ with a > 0 does not belong to the same universality class as the KSE ([^. 
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